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We examine, using molecular dynamics simulation, the structure and thermodynamics of the (100) 
and (111) disordered face-centered cubic (FCC) crystal/melt interfaces for a binary hard-sphere 
system. This study is an extension of our previous work, [Phys. Rev. E 54, R5905 (1996)], in 
which preliminary data for the (100) interface were reported. Density and diffusion profiles on 
both fine- and course-grained scales are calculated and analyzed leading to the conclusion that 
equilibrium interfacial segregation is minimal in this system. 



I. INTRODUCTION 

In order to fully understand such important phenom- 
ena as near-equihbrium crystal growth, homogeneous nu- 
cleation and interfacial solute segregation, a detailed mi- 
croscopic description of the crystal/melt interface is nec- 
essary . Since direct experimental data is scarce, due 
to the extreme difficulty of constructing and interpreting 
experiments on such systems, computer simulation has 
become an important tool, not only in its usual role of 
aidingin the development and evaluation of interface the- 
ories [Q-D , but also in determining the basic microscopic 
phenomenology of crystal/melt interfaces. Previous sim- 
ulation studies of such interfaces have focused entirely on 
single-component systems (for a review of recent simula- 
tion studies, see Ref. [||); however, since most technolog- 
ically important materials are not pure substances but 
mixtures (such as alloys), it is important that such stud- 
ies be extended to multi-component interfacial systems. 
In addition to the usual issues of interfacial width, in- 
terfacial free energy, and transport within the interface, 
multi-component systems allow for the study of interfa- 
cial segregation. 

As in single-component systems, the crystal/melt in- 
terface of a multi-component system is characterized by 
measuring the change in the various structural, ther- 
mal and dynamical properties of interest as one traverses 
the interface from one bulk phase into the other. For 
planar interfaces (the type studied here), the z axis is 
usually taken as the direction perpendicular to the in- 
terfacial plane and quantities are averaged over x and 
y and presented as functions of z. Examples include 
the density profiles of the various components (labeled 
by i), Pi{z) = {pi{r))xy and the diffusion constant pro- 
files, Di{z). The thermodynamic quantity of most inter- 
est is the solid- liquid interfacial free energy, 7si, which 
is defined as the work required to form one unit area 
of surface - this quantity is extremely difficult to calcu- 
late via simulation - even in the single-component case. 
The only reliable calculation is the calculation of 7si for 
a single-component Lennard- Jones system by Broughton 
and Gilmer |^]. 



Recently we have reported initial simulation re- 
sults on the crystal/melt interface of a two-component 
hard-sphere mixture. Specifically, the interface between 
the (100) face of a disordered face-centered cubic (FCC) 
crystal and the coexisting fluid was studied. In this work 
we extend this calculation to the (111) face and revisit the 
(100) interface calculations using more detailed analysis. 
These simulations represent the first such simulations on 
the crystal/melt interface of a multi-component system. 

The reasons for choosing the hard-sphere system for 
this initial study are two-fold: First, it is now well estab- 
lished that the structure and freezing behavior of dense, 
simple fluids, is determined, for the most part, by pack- 
ing considerations determined by the repulsive part of the 
interaction potential. The effect of the attractive forces 
can generally be accounted for by treating it as a per- 
turbation to the repulsive part of the potential, which is 
often approximated by a hard-sphere with some effective 
diameter Thus, the hard-sphere system is a useful 
reference from which to begin studies of more realistic 
systems. Studying the hard-sphere system also allows 
one to directly probe the role of packing (which is purely 
entropic) in determining the interfacial phenomenology. 
Second, the relative simplicity of the hard-sphere system 
lends itself well to theoretical study ~ the vast majority of 
density- functional calculations on crystal/melt interfaces 
involve hard-sphere systems. The disordered-FCC/fluid 
interfaces were chosen for initial study, because, in or- 
der to begin an interface simulation, the phase coexis- 
tence conditions must be very accurately known, and the 
disordered-FCC/fluid region of the binary hard-sphere 
phase diagram (which occurs in the region of the phase 
diagram where the difference between the diameters of 
the different components is not too great) has been well 
characterized pO| ]. 

In the next section we describe the binary hard-sphere 
system and its phase behavior. The results of the disor- 
dered FCC/fluid interfacial simulations reported for the 
(100) and (111) interfaces are discussed in Section 3. In 
Section 4, we conclude. 
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II. THE BINARY HARD-SPHERE SYSTEM 

We consider a two-component system consisting of 
hard spheres of differing sizes. The interaction potential 
for such a system can be written 

where r is the distance between the centers of two spheres 
and i and j £ {1,2} index the two types of spheres, 
which are distinguished by their different diameters, cri 
and {72, while their masses are assumed to be identical. 
We also assume that the spheres are additive; i.e., tJn = 
fi, o'22 = o'2, and 0-12 = cr2i = (cti -I-(T2)/2. The following 
definitions and conventions are adopted for this study. 
First, it is assumed, without loss of generality, that <72 > 
(Ti, and the hard-sphere diameter ratio a is defined as 
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a = oxja-i, (0 < a < 1) . 



(2) 



Second, if there are A'^i hard spheres with diameter a\ 
and N2 with diameter 02 in the volume V ^ then 



A^i + A^2 

P = y = PI+P2 



(3) 



is the total number density, and pi's represent the respec- 
tive number densities for individual species. Third, the 
concentrations (mole fractions) of individual components 
are given by 



Pt/P, i = l,2 



(4) 



Since xi + X2 — 1, a single variable x is usually used, 
such that X2 = X and x^ = \ — x. 

Also, the total packing fraction for the mixture 77 is 



77 = r/i + r\2 



(5) 



where r\i = '^a\ pi, « = 1,2. The diameter ratio a to- 
gether with any pair of independent parameters from 
those defined above can be used to completely specify 
the fluid state of a binary hard-sphere system. The unit 
of length for the binary system is taken to be equal to 
the diameter of a larger sphere 02- 

Depending on the value of the diameter ratio a, dif- 
ferent crystal structures may have the lowest free en- 
ergy |l^Jl^. In particular, for a above about 0.85, the 
substitutionally disordered FCC crystal, in which spheres 
of different diameters are distributed randomly over the 
sites of an FCC lattice, is the stable structure at freezing. 
For a below this value, a variety of structures may ex- 
ist including ordered solid such as NaCl and CsCl struc- 
tures iQ. For an interface simulation study, it is neces- 
sary that the phase coexistence conditions be accurately 
known. Therefore, we have chosen for this study the 
disordered-FCC/melt system, since the coexistence con- 
ditions for the equilibrium disordered crystal/melt in- 
terface can be found in the study of Kranendonk and 



FIG. 1. Phase diagram of a binary hard-sphere system with 
a = 0.9 . The data are taken from simulation by Kranendonk 
and Frenkel . The triangles correspond to the fluid phase, 
and circles to the crystal phase. The solid triangle and cir- 
cle represent the coexistence state selected for the interface 
simulation in the present work. 



Frenkel [|lO|, who have calculated the crystal/melt phase 
diagrams for several values of the diameter ratio in the 
range 0.85 < a < 1. For the present study, the diameter 
ratio a = 0.9 is selected. 

For a binary system with a given diameter ratio a, the 
coexistence state is specified by two number densities p'^ 
and p\ as well as two concentrations x"^ and x^ (or, alter- 
natively, by four number densities of individual species, 
p^ and pI, where i = 1,2). The pressure-concentration 
phase diagram determined in the simulation by Kranen- 
donk and Frenkel ||l^ for a = 0.9 is shown in Fig. 0. At 
a given pressure, the fluid and crystal phases have dif- 
ferent concentrations at coexistence, represented in the 
plot by triangles and circles, respectively. To maximize 
the deviation from single-component behavior, a point on 
the a = 0.9 phase diagram was chosen for the interface 
simulation, where the concentration difference between 
crystal and melt phases is the largest. This point occurs 
at a pressure of 



P = 14.7A:bT/ct| 



and concentrations of 



= 0.71 and 



= 0.54 



(6) 



(7) 



for the crystal and fluid phase, respectively, and is shown 
in Fig. |l| by the solid triangle and circle. We have found 
the total packing fractions (number densities) in each two 
phases at 



= 0.552 
= 0.502 



(pV| = 1.144), 
(p'cr| = 1.096) 



(8) 
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and have used these values for the initiahzation of the 
binary interface systems. We have run simulations for 
several trial systems and have found that the systems are 
stable and the bulk crystal remains on average stress- free. 
Therefore, the above parameters have been used for all 
the simulation runs from which average properties of the 
interface have been computed. 

III. INTERFACE CONSTRUCTION AND 
SIMULATION RESULTS 

To create the initial bulk systems that are placed to- 
gether to form the interface, the two sphere types are 
distributed randomly according to the crystal and fluid 
coexistence concentrations. The concentration in each 
crystal layer is maintained fixed by randomly distribut- 
ing the spheres of different types on a layer by layer ba- 
sis, thereby removing layer-to-layer concentration fluctu- 
ations due to finite system size. This constraint is not 
expected to affect the results in any significant way be- 
sides removing the fluctuations that cannot be averaged 
over during the simulation run due to practically no dif- 
fusion in the bulk crystal. The random distribution of 
particle types is justified by the conclusion of Kranen- 
donk and Frenkel that above about x = 0.6 there is little 
or no local substitutional ordering |]l3| . 

Since the hard-sphere system evolves on a collision-by- 
collision basis, the natural unit of time for hard spheres is 
the mean collision time Tc, i.e. the average time between 
collisions suffered by a given particle. On the other hand, 
the duration of a simulation run is most conveniently 
measured in terms of the total number of collisions. In 
the present study, in order to have better correspondence 
between the simulation time and the system evolution 
time, we measure simulation time in units of the number 
of collisions per particle (cpp), defined as twice the ratio 
of the total number of collisions to the number of particles 
in the system, so that 1 cpp « Tc . 

We have prepared 10 systems for each of (100) and 
(111) crystal orientations and have run the simula- 
tions for 20 000 cpp with the interfacial diagnostics be- 
ing recorded every 200 cpp. The (100) systems contain 
11616 spheres and have dimensions = Ly — 16.70(72, 
and Lz — 37.09(72, while the size of the (111) systems 
is 11340 spheres and Lx = 16.09(72, Ly — 16.72(72, and 
Lz = 37.60(72. (We have also done simulations on these 
systems with smaller numbers of spheres ~ 3000 and 
6000 - and different cross-sectional shapes. The results 
of these simulations are within statistical error, quan- 
titatively identical to those presented here. The larger 
samples have been chosen as they give better statistics, 
have much shorter equilibration times than the smaller 
samples, and exhibit smaller interfacial fluctuations.) 

The concentration fluctuations in the fluid and the in- 
terfacial regions have been found to be much larger and 
more persistent than the density fluctuations. Also, due 
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FIG. 2. Fine-scale density profiles for the (100) binary mix- 
ture interface. 

to a finite system size, even though the total momen- 
tum with respect to the simulation cell is set to zero, 
drift of the interface positions is observed. In order to 
avoid broadening of the interfacial profiles caused by the 
drift, we have selected for the final averages 12 segments 
of 2000 cpp in duration from each crystal orientation, 
such that the drift does not exceed half the distance be- 
tween crystal layers (for more details on the methods 
of interfacial construction and equilibration used here, 
as well as more detailed definitions of measured (com- 
puted) profiles, see our recent work on single-component 
systems Q). 

The fine-scale profiles for the two components pi{z), 
P2{z) and for the total density p{z) = pi{z) + P2(z) are 
shown in Figs. ^ and || for the (100) anci (111) crystal 
orientations, respectively. For the total density profiles 
the 10-90 widths of the height of the density peaks equal 
to 5.3 and 5.6(72 for the (100) and (111) interfaces, re- 
spectively. The density oscillations in p2{z) and p{z) 
dampen monotonically, while pi{z) exhibits a peculiar 
non-monotonic peak-height envelope, a phenomenon that 
has not been seen in any of the single-component system 
studies. 

The non-monotonic behavior of the fine-scale density, 
pi , can be explained by examining the coarse-grained (fil- 
tered) density profiles pi{z), P2{z), and p{z), computed 
using Finite Impulse Response (FIR) filter |l3| . The use 
of such filters for coarse-graining the density profiles is 



3 



I J 


1 - 


1 1 




1 1 


1 / 




u 


i 


(111) 


1 j 


11 




I 


1 , 


1 1 


, 1 


I- 


WW 










1 






1 


I 




, , , ^ 



-8 -6 -4 -2 2 4 6 8 

2;, 0"2 

FIG. 3. Fine-scale density profiles for the (111) binary mix- 
ture interface. 

necessary when the peak-to-peak spacing of a profile is 
not constant through the interface - in such a case the 
use of uniform bins to perform coarse graining can lead 
to misleading results. (For a detailed description of the 
use of such filters in analyzing density profiles, see our 
recent article on the single-component hard-sphere in- 
terface 10.) The filtered profiles are shown in Figs. ^ 
and|| for the (100) and (111) interfaces, respectively. In 
both cases the individual species densities change have 
10-90 widths of about 2(72, whereas the corresponding 
width for the total density is about 4a2. This seems 
strange in that the total density is the sum of the in- 
dividual densities, but since the difference between the 
total densities on either side of the interface is an or- 
der of magnitude smaller than that for the individual 
densities, very small changes in pi or p2 can contribute 
significantly to the 10-90 width of the total density, while 
remaining unimportant in determining the width of the 
individual densities. (Note that, the width of the total 
density is somewhat larger than the 3.2-3.3(7 found in the 
single-component case |14| ). The rapid change in concen- 
tration over about 2(72, combined with the fact that the 
density oscillations on the fluid side of the interface in 
the fine-scale profiles persist 2-3(72 after the concentra- 
tions have relaxed, leads to the above mentioned non- 
monotonicity. Since the number density of the smaller 
spheres in the fluid [p{ = (1 - x^)p^ = 0.504(7^^] is 
larger than the corresponding density in the crystal re- 
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FIG. 4. Filtered density profiles for the (100) binary mix- 
ture interface. 
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FIG. 5. Filtered density profiles for the (111) binary mix- 
ture interface. 

gion [pI = (1 — x'^)p'^ = 0.332 (72^'^], the ordering of the 
fluid in the presence of the interface occurs in a region 
with higher average density than that in the bulk crys- 
tal, resulting in the higher proflle peaks in the interfacial 
region. 

Analysis of the coarse-grained densities also leads to 
the conclusion that there is no statistically signiflcant 
equilibrium interfacial segregation (Gibbs adsorption) in 
either of our simulated interfacial orientations. Such seg- 
regation is quantifled by Fi/A, where Fi is the excess 
number of type 1 (small) particles (here the smaller par- 
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tide is taken to be the solute), defined using a Gibbs di- 
viding surface |^,|^ for which the excess number of type 2 
(large) is zero, and A is the area of the interface. The sig- 
nificance of this result should be taken, however, in light 
of the issues of chemical equilibrium discussed below. 

Information about changes in the inter-layer spacing is 
obtained by measuring the layer separation, defined as 



Az,, 



(9) 



where Zi is the center of mass of layer i determined from 
the fine-scale density profile between the adjoining den- 
sity minima. We calculate Az^ for the total density pro- 
file as well as for the density profiles of individual com- 
ponents. The results are shown in Fig. |^ with diamonds, 
squares and circles representing layer separation in pi{z)^ 
P2{z) and p(z), respectively. As in the single-component 
case [0, the layer separation shows large layer expan- 
sion for the (100) interface and very little expansion for 
the (111) interface. In addition, we see significantly dif- 
ferent behavior of the individual density profiles. For the 
(100) orientation the increase in the layer separation of 
the pi{z) profile is delayed by about two crystal layers, 
compared with that of P2{z)- For the (111) orientation 
the pi{z) profile exhibits almost no change in the inter- 
layer spacing. Evidently, at the onset of the disorder 
in the interfacial region, the spheres of type 2, having 
larger diameter, are repelled farther away from the or- 
dered crystal layers. 

The self-diffusion constant profiles are computed sepa- 
rately for the two particle types according to the formula 
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FIG. 7. Diffusion profiles for smaller (circles) and larger 
(asterisks) spheres in the binary (100) and (111) crystal/melt 
interfaces. The error bars represent twice the standard devi- 
ation of the mean value calculated from the 12 samples. 
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where Ni{z) is the number of spheres of type i between 
z — Sz/2 and z + Sz/2 at time t — to, where Sz is the 
crystal layer spacing, and the brackets represent the av- 
erage over time origins to. The sphere displacement was 
monitored on a uniform coarse scale over time imax — = 
5.5 (mcrl/fceT)^/^. The averaging was done over 50 time 
origins for each of the 12 selected intervals. The average 
diffusion constant profiles for the (100) and (111) inter- 
faces are shown in Fig. |^. As could have been anticipated, 
the self-diffusion coefficient in the bulk fluid is larger for 
the smaller spheres at 0.019 (A;bT(t|/to)^/^, compared to 
that for the larger spheres at 0.016 (fcBTcrl/m)^/^. The 
diffusion profiles change monotonicly across the inter- 
face with the 10-90 widths of about 3.3 cr2 for both (100) 
and (111) interfaces, which is intermediate between the 
widths for the interfacial profiles and that for the total 
density. Note that, since the midpoint of the diffusion 
profile is shifted by about 1-2(72 to the liquid side of the 
midpoints of the density profiles (that is, the bulk of the 
increase in the diffusion constant occurs after the densi- 
ties have relaxed to nearly their liquid values), the actual 
width of the interfacial region is larger than either trans- 
port or structural properties would indicate alone. 

At this point it is useful to comment on the question 
of chemical equilibrium in our systems. When we create 
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the interface by placing crystal and fluid blocks next to 
each other and then allowing the system to evolve, we ex- 
pect that after some time the system will stabilize itself 
and the interfaces will be in equilibrium. In our sim- 
ulations we see that after the initial equilibration, both 
temperature (see Ref. j|] ) and transverse pressure profiles 
exhibit no significant deviation in the interfacial region 
from their average values. This is an indication that the 
interface is in thermal and dynamic equilibrium with the 
surrounding bulk phases. On the other hand, the concen- 
tration equilibrium of individual species at the interface 
cannot be assumed with the same certainty. One admit- 
ted approximation that we have made in the construction 
of our interface is the use of a randomly substituted solid 
mixture. In their studies on hard-sphere binary mixtures, 
Kranendonk and Frenkel |13| report no significant local 
substitutional order for solid solutions above 60% large 
spheres. Since our solid has a large-sphere mole fraction 
of 0.71, our assumption of random substitution in the 
bulk solid should be valid. However, the work of Kranen- 
donk and Frenkel applies strictly to the bulk system and 
it is possible that there is some equilibrium local substi- 
tutional order that should be present on the solid side of 
our interface that we cannot see due to the very long re- 
laxation times for concentration fluctuations. The closer 
we are to the crystal side of the interface, the less certain 
can we be that a particular configuration of the small and 
large spheres correctly represents the equilibrium concen- 
tration profile. The reason is that, unlike temperature 
and pressure fluctuations which are transported through 
the system via collisions, the concentration fluctuations 
are introduced when particles themselves drift through 
the system. This obviously requires moderate values of 
the diffusion coefficient which cannot be achieved in the 
crystal. (In order to achieve true concentration equilib- 
rium in a binary interfacial system one would need to sim- 
ulate composition fiuctuations consistent with the chem- 
ical potential balance in both crystal and fluid phases 
and across the interface. This can be done, for example, 
by introducing Monte-Carlo moves into the equilibration 
process that allow small spheres to become large ones 
and vice versa, with probabilities that produce correct 
chemical potential profiles. We have not, however, found 
a consistent way of doing this in practice that both pre- 
serves the stability of the interface and gives good statis- 
tics.) This being said, the effect of the uncertainty in the 
degree of chemical equilibrium achieved should have little 
effect on most of the phenomena mentioned above, such 
as non-monotonicity of the pi profile and the anomalous 
behavior of the lattice spacing in both interfaces, since 
these effects are primarily due to behavior on the liquid 
side of the interface where the diffusion constant is large 
enough to ensure relaxation in concentration. The largest 
effect will be on the width of the concentration profiles 
on the solid side of the interface - the value given above 
should be taken as an upper bound - and the precise de- 
gree of equilibrium solute segregation (adsorption) on the 
solid side. The high stability of our interfaces leads us to 



speculate that these effects will be small, but, because of 
the problems outlined above, they cannot be discounted. 



IV. SUMMARY AND CONCLUSIONS 

We have presented detailed molecular-dynamics sim- 
ulations for the (100) and (111) (disordered) FCC crys- 
tal/melt interfaces for a binary system of hard spheres. 
The ratio of the small hard-sphere diameters and that of 
the large spheres was chosen to be 0.9. This study ex- 
tends our earlierpreliminary work on the (100) interface 
for this system g . The principle results of this study are 
summarized as follows: 

1. The fine-scale density profiles for the smaller par- 
ticle, pi{z)^ in contrast to the single-component 
case, exhibit a pronounced non-monotonic enve- 
lope. This behavior is not seen in either the to- 
tal or large-particle density profiles. Analysis of 
the coarse-grained density profiles shows that this 
phenomenon is not due to any appreciable adsorp- 
tion of the smaller particle at the interface, but is 
entirely due to the fact that increase in the small 
particle concentration occurs over a shorter length 
scale than the decay of the density ocsillations in 
the liquid. 

2. As in the single-component case |^ we see an in- 
crease in the spacing between the (100) density- 
profile peaks as the interface is traversed from crys- 
tal to fluid. A much less pronounced effect is seen in 
the (111) interface where the large-particle density- 
peak spacing stays mostly constant except for a 
maximum well on the liquid side; in contrast, the 
spacing for the smaller particles is constant through 
the interface. 

3. The widths of the coarse-grained concentration pro- 
flies (calculated with a Finite Impulse Response 
(FIR) filter) are considerably smaller that those 
for the total densities (about 2(72 versus 4(T2) and 
no significant equilibrium interfacial segregation is 
seen. As in our earlier single-component study , 
the use of FIR filters to determine coarse-grained 
density profiles is necessary to avoid artifacts of the 
binning process when the peak separation is not 
constant. 

4. The diffusion profiles and coarse-grained density 
are shifted by 1 to 2ct2 relative to one another. 
Significant diffusion begins only after the density 
has relaxed to nearly the bulk liquid value. There- 
fore, the interfacial region is wider than either of 
these profiles would indicate separately, and we can 
identify two separate sub-regions as the interface 
is traversed from solid to fluid, in which density 
relaxation or changes in transport properties are 
dominant, respectively. 
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